#### 准备好环境 ####

#### 1 定位数据模块和变量,获取数据源 ####
# 2023年2月17日20:26:03 明天开始正式学习

# 1.1 获取demo_g文件 这个包加载顺序有次序 暂时还不知道先加载谁,看视频没看出问题 需要加载两次
library(haven) #用于读取xpt
library(plyr) #用于数据处理包
library(dplyr) #用于数据处理包
library(arsenal) ##用于tableby
library(survey) #用于加权分析
library(grid)
library(Matrix)

setwd("G:/BaiduNetdiskDownload")
#### 1.1 DEMO-人口学数据提取 #####

demo.g <- read_xpt("NHANES/2011-2012/Demographics/demo_g.xpt")
demo.h <- read_xpt("NHANES/2013-2014/Demographics/demo_h.xpt")

#colnames(demo.g)
#colnames(demo.h)
### 1.2 提取研究所需要的变量 
#tableby是用于快速查看表数据用的 可以用因子处理类似性别的数据 快速查看
tab1 <- tableby( ~ RIDAGEYR + factor(RIAGENDR) + (RIDRETH3) + DMDEDUC2 + INDFMPIR,#需要显示的列
                #strata = factor(RIAGENDR),
                data = demo.g)
 summary(tab1,text = T)
# 部分数据出现了n-miss 是指的数据丢失 导致 元数据补全 可以通过查看原始数据得到

# demo.g_colnames<-colnames(demo.g)
# colnames(demo.g)
# colnames(demo.h)
# aaa<- as.array(demo.g_colnames)
# class(aaa)
# print(aaa)
# list(aaa)
# ?remove

# demo.g_colnames["INDFMPIR"]
# print(demo.g)
# |                                  | Overall (N=9756) |
#   |:---------------------------------|:----------------:|
#   |Age in years at screening         |                  |
#   |-  Mean (SD)                      | 31.403 (24.579)  |
#   |-  Range                          |  0.000 - 80.000  |
#   |factor(RIAGENDR)                  |                  |
#   |-  1                              |   4856 (49.8%)   |
#   |-  2                              |   4900 (50.2%)   |
#   |Race/Hispanic origin w/ NH Asian  |                  |
#   |-  Mean (SD)                      |  3.440 (1.602)   |
#   |-  Range                          |  1.000 - 7.000   |
#   |Education level - Adults 20+      |                  |
#   |-  N-Miss                         |       4196       |
#   |-  Mean (SD)                      |  3.467 (1.283)   |
#   |-  Range                          |  1.000 - 9.000   |
#   |Ratio of family income to poverty |                  |
#   |-  N-Miss                         |       840        |
#   |-  Mean (SD)                      |  2.206 (1.638)   |
#   |-  Range                          |  0.000 - 5.000   |
#年龄- RIDAGEYR 性别-> RIAGENDR  种族-> RIDRETH3 教育程度-> DMDEDUC2 贫困程度->INDFMPIR 
#两个数据合并显示 根据行进行组合一起 指的是上下堆叠合并一个大表 有点类似 123 + 456 行累加起来!

demo.data.file <- dplyr::bind_rows(list(demo.g,demo.h))
#View(demo.g)
#View(demo.h)
#View(demo.data.file)

demo.data <- demo.data.file[,c('SEQN', 'RIDAGEYR', 'RIAGENDR', 'RIDRETH3', 'DMDEDUC2', 'INDFMPIR')]
#View(demo.data)#查看数据
# 分类变量用factor
#strata  分组用的吗? 这块我看视频没起作用 通过给出不同的字段进行分组了数据 确实很厉害这块 如果是java 需要搞很麻烦的写法
#data 是数据来源 tableby 是查看部分变量分布情况 加上strata 是进行根据分组显示 factor 是进行分类变量 使用 类似 性别用的  

# tab1 <- tableby( ~ RIDAGEYR + factor(RIAGENDR) + factor(RIDRETH3) + factor(DMDEDUC2) + INDFMPIR, 
#                  data=demo.g)  
# summary(tab1, text=TRUE)

##### 1.2 SMQ-吸烟数据提取 #####
### 1.1 提取 Component文件
smq_g <- read_xpt("NHANES/2011-2012/Questionnaire/smq_g.xpt")
smq_h <- read_xpt("NHANES/2013-2014/Questionnaire/smq_h.xpt")

### 2. 提取研究所需要的变量
# 是否吸烟至少100支-SMQ020; 现在是否吸烟-SMQ040; 多久前开始戒烟-SMQ050Q;
# 多久前开始戒烟的时间单位（天、周、月、年）-SMQ050U;

smq.data.file <- dplyr::bind_rows(list(smq_g,smq_h))
#View(smq.data.file)
colnames(smq.data.file)
smq.data <- smq.data.file[,c('SEQN', 'SMQ020', 'SMQ040', 'SMQ050Q', 'SMQ050U')]
# View(smq.data)
##### 1.3 ALQ-饮酒数据提取 ######

alq_g <- read_xpt("NHANES/2011-2012/Questionnaire/alq_g.xpt")
alq_h <- read_xpt("NHANES/2013-2014/Questionnaire/alq_h.xpt")

### 2.提取研究所需要的变量(木)
# 每年至少喝12杯酒-ALQ101
# 一生中至少喝过12次酒-ALQ110；
# 过去12个月多久喝一次酒-ALQ120Q
# 过去12个月饮酒频率的单位（周，月，年）-ALQ120U;
alq.data.file <- dplyr::bind_rows(list(alq_g, alq_h))
# colnames(alq.data.file)
alq.data <- alq.data.file[,c('SEQN', 'ALQ101', 'ALQ110', 'ALQ120Q', 'ALQ120U')]
#View(alq.data)
# tab1 <- tableby(~factor(ALQ101) + factor(ALQ110) + ALQ120Q + factor(ALQ120U),
#                data = alq.data)
# summary(tab1, text=TRUE)

##### 1.4 BMX-BMI、腰围数据提取 #####
### 1.提取 Component文件
#注意是 Examination 的类别
bmx_g <- read_xpt("NHANES/2011-2012/Examination/bmx_g.xpt")
bmx_h <- read_xpt("NHANES/2013-2014/Examination/bmx_h.xpt")

bmx.data.file <- dplyr::bind_rows(list(bmx_g, bmx_h))
bmx.data <- bmx.data.file[,c('SEQN', 'BMXBMI', 'BMXWAIST')]
#View(bmx.data)


##### 1.5 饮食-L&Z 以及总卡路里数据提取(第1天 & 第2天) #####

# 第1天
dr1tot.g <- read_xpt('NHANES/2011-2012/Dietary/dr1tot_g.xpt')
dr1tot.h <- read_xpt('NHANES/2013-2014/Dietary/dr1tot_h.xpt')
# View(dr1tot.g)

dr1tot.data.file <- dplyr::bind_rows(list(dr1tot.g, dr1tot.h))
dr1tot.data <- dr1tot.data.file[,c('SEQN', 'DR1TLZ', 'DR1TKCAL')]
# View(dr1tot.data)

# 第2天
dr2tot.g <- read_xpt('NHANES/2011-2012/Dietary/dr2tot_g.xpt')
dr2tot.h <- read_xpt('NHANES/2013-2014/Dietary/dr2tot_h.xpt')
dr2tot.data.file <- dplyr::bind_rows(list(dr2tot.g, dr2tot.h))
dr2tot.data <- dr2tot.data.file[,c('SEQN', 'DR2TLZ', 'DR2TKCAL')]
# View(dr2tot.data)
#是把两个表安列合并 list 是上下 merge 是左右合并 根据唯一标识符 seqn 
dr.data <- merge(dr2tot.data,dr1tot.data)

# View(dr.data)

##### 1.6 膳食补充剂-L&Z数据提取(第1天 & 第2天) ######
# 第1天
ds1tot.g <- read_xpt('NHANES/2011-2012/Dietary/ds1tot_g.xpt')
ds1tot.h <- read_xpt('NHANES/2013-2014/Dietary/ds1tot_h.xpt')
# View(ds1tot.g)
ds1tot.data.file <- dplyr::bind_rows(list(ds1tot.g, ds1tot.h))
ds1tot.data <- ds1tot.data.file[,c('SEQN', 'DS1TLZ')]


# 第2天
ds2tot.g <- read_xpt('NHANES/2011-2012/Dietary/ds2tot_g.xpt')
ds2tot.h <- read_xpt('NHANES/2013-2014/Dietary/ds2tot_h.xpt')
ds2tot.data.file <- dplyr::bind_rows(list(ds2tot.g, ds2tot.h))
ds2tot.data <- ds2tot.data.file[,c('SEQN', 'DS2TLZ')]

# 合并第1天&第2天的变量
ds.data <- merge(ds1tot.data, ds2tot.data)
# View(ds.data)

##### 1.7 CFQ-认知数据提取 #####
cfq.g <- read_xpt('NHANES/2011-2012/Questionnaire/cfq_g.xpt')
cfq.h <- read_xpt('NHANES/2013-2014/Questionnaire/cfq_h.xpt')

### 2.提取研究所需要的变量
#CERAD 词汇第1次即刻记忆测验	CFDCST1
#CERAD 词汇第2次即刻记忆测验	CFDCST2
#CERAD 词汇第3次即刻记忆测验	CFDCST3
#CERAD 3次词汇即刻记忆测验总分	CERAD_1+CERAD_2+CERAD_3
#CERAD 延迟回忆评分	CFDCSR
#语言流畅性评分Animal Fluency	CFDAST
#数字符号替代测试(DSST)	CFDDS

cfq.data.file <- dplyr::bind_rows(list(cfq.g, cfq.h))
cfq.data <- cfq.data.file[,c('SEQN', 'CFDCST1', 'CFDCST2', 'CFDCST3',
                             'CFDCSR', 'CFDAST', 'CFDDS')]
#View(cfq.data)



#### 2 提取分析相关变量（权重等，暂时为复现paper结果而提取） ####
##### 2.1 权重变量 ##### 
dr1tot.data.add.weight <- dr1tot.data.file[,c('SEQN', 'WTDRD1', 'WTDR2D', 'DR1TLZ', 'DR1TKCAL')]
dr2tot.data.add.weight <- dr2tot.data.file[,c('SEQN', 'WTDRD1', 'WTDR2D', 'DR2TLZ', 'DR2TKCAL')]
dr.data.add.weight <- merge(dr1tot.data.add.weight, dr2tot.data.add.weight)
#View(dr.data.add.weight)


weight.data <- dr1tot.data.file[,c('SEQN', 'WTDRD1')]
weight.data$WTDRD1 <- weight.data$WTDRD1/2 #权重的计算方式，需要按 Cycle 进行平均
#View(weight.data)


##### 2.2 复杂抽样的其他变量-DEMO #####
survey.design.data <- demo.data.file[,c('SEQN', 'SDMVPSU', 'SDMVSTRA')]
#View(survey.design.data)

#### 3 合并上述所有步骤 ####
##### 3.1 合并步骤1 & 2中提取的数据 #####
#demo.data, smq.data, alq.data, bmx.data, dr.data, ds.data, cfq.data, weight.data
# plyr::join_all 是合并各个分数据到一个大表 通过特定的主键 左连接方式
output <- plyr::join_all(list(demo.data, smq.data, alq.data, bmx.data,
                              dr.data, ds.data, cfq.data, weight.data, survey.design.data),
                         by='SEQN', type='full')
# View(output)
# dim(output)
#筛选单元格根据值进行拆分

##### 3.2 针对合并后的数据进行质控，确定人数能对得上 #####
data.age.60 <- subset.data.frame(output, RIDAGEYR >= 60 )

# dim(data.age.60) #[1] 3632    

# 2. 确认下权重 挑选不是na 的也就是没有有效数据的过滤掉
which(!is.na(data.age.60$WTDRD1))
weight.non.na.index <- which(!is.na(data.age.60$WTDRD1))
# 计算权重-59659564 gcr 表格权重 数值可以对上 这句话是表格里面WTDRD1 权重 不是na的行数拿出来取和
sum(data.age.60$WTDRD1[weight.non.na.index]) 

#### 4.根据文章的筛选策略进行筛选 ####
# 原始筛选策略：The analyses presented here include a pooled sample of all individuals aged
# 60 years or older in each of the two survey cycles, with completed cognitive performance
# test results, dietary L and Z intake information, and information on important
# con-founders(age, sex, race/ethnicity, smoking, educational attainment).
paper.data <- subset.data.frame(output, RIDAGEYR >= 60 &
                                  (!is.na(CFDCSR)) & #CERAD 延迟回忆评分
                                  (!is.na(CFDAST)) & #语言流畅性评分Animal Fluency
                                  (!is.na(CFDDS)) & #数字符号替代测试(DSST)
                                  (!is.na(DR1TLZ))& #饮食-叶黄素 & 玉米黄质摄取
                                  (!is.na(DR1TKCAL))& #能量摄入总量-Energy (kcal)
                                  (!is.na(RIDAGEYR)) & #年龄
                                  (!is.na(RIAGENDR)) & #性别
                                  (!is.na(RIDRETH3)) & #种族
                                  (!is.na(DMDEDUC2)) & #教育程度
                                  (!is.na(SMQ020))) #是否吸烟至少100支)
 
# 上面的取值会影响下面每个维度的统计百分比,一方面是数据可能有变化 另一方面是加权可能不一致导致 目前是total 这个维度是一致的 纯粹大于60 以上的人
# paper.data <- subset.data.frame(output, RIDAGEYR >= 60 
# )

# 这个值与n= 2796 对不上可能是数据右边或者统计其他维度没给说 目前这块是根据木统计得出
dim(paper.data)


#View(paper.data)
#数据可能对不上是因为数据可能会变化 原来统计错误等 会修正错误 2023年2月20日11:06:57
paper.exclude.data <- subset.data.frame(output, RIDAGEYR >= 60 &(
  (is.na(CFDCSR)) | # CERAD 延迟回忆评分
    (is.na(CFDAST)) |# 语言流畅性评分Animal Fluency
    (is.na(CFDDS)) | # 数字符号替代测试(DSST)
    (is.na(DR1TLZ))| #饮食-叶黄素 & 玉米黄质摄取
    (is.na(DR1TKCAL))| #能力摄入总量-Energy (kcal)
    (is.na(RIDAGEYR)) | # 年龄
    (is.na(RIAGENDR)) | # 性别
    (is.na(RIDRETH3)) | # 种族
    (is.na(DMDEDUC2)) | # 教育程度
    (is.na(SMQ020))) # 是否吸烟至少100支
)
# dim(paper.exclude.data)
# colnames(paper.exclude.data)
#$Type 是打上标签 相当于在表头上面加一个开头
paper.exclude.data$Type <- 'exclude'
paper.data$Type <- 'include'
compare.data <- dplyr::bind_rows(paper.data, paper.exclude.data)
# paper.data$Type <- 'NHANES participants included in this analysis(n=2796)'
#下面这个是我新加的
# data.age.60$Type <- 'All NHANES participants aged ≥60 years(n=3632)'
# dim(data.age.60)
# compare.data <- dplyr::bind_rows(paper.data,paper.exclude.data,data.age.60)
# compare.data <- dplyr::bind_rows(paper.data,data.age.60)
# 后续统计可以吧新的compare.data 赋值给paper.data
# paper.data <- compare.data


#tabexcludeinclude <- tableby(Type ~ factor(DMDEDUC2), data = paper.data) # 复现论文中的排除人群和纳入人群在人口学变量上的差异


##### 4.1.2 进行实现论文表格中数据表的复现开始复现表格数据中间 ##### 
## 进行复现数据展示2023年2月21日16:09:29
paper.data$Sex <- ifelse(paper.data$RIAGENDR==1,'male','female')
# length(paper.data)
# summary(paper.data)


# compare.data$Sex <- ifelse(compare.data$RIAGENDR==1,'male','female')
#  length(compare.data)
# summary(compare.data,text=T)
# 
# compare.datatableSexView<- tableby(Type ~ factor(Sex),data = compare.data,weights = WTDRD1,subset = WTDRD1 >0)
# summary(compare.datatableSexView,text=T)

tableSexView<- tableby(Type ~ factor(Sex),data = paper.data,weights = WTDRD1,subset = WTDRD1 >0)
summary(tableSexView,text=T)
# |            | include (N=54660094) | Total (N=54660094) |
#   |:-----------|:--------------------:|:------------------:|
#   |factor(Sex) |                      |                    |
#   |-  female   |   29455555 (53.9%)   |  29455555 (53.9%)  |
#   |-  male     |   25204539 (46.1%)   |  25204539 (46.1%)  |

##### 4.1.3 年龄-age.group #####  
paper.data$age.group <- ifelse(paper.data$RIDAGEYR >= 60 & paper.data$RIDAGEYR < 69, '60-69 years',
                               ifelse(paper.data$RIDAGEYR >=70 & paper.data$RIDAGEYR <79, '70-79 years',
                                      '80+ years'))
table(paper.data$age.group)

tab1 <- tableby(Type ~ factor(age.group),
                data = paper.data, weights=WTDRD1,
                subset = WTDRD1>0) 

summary(tab1, text=TRUE) 

# 
# |                  | include (N=54660094) | Total (N=54660094) |
#   |:-----------------|:--------------------:|:------------------:|
#   |factor(age.group) |                      |                    |
#   |-  60-69 years    |   28624430 (52.4%)   |  28624430 (52.4%)  |
#   |-  70-79 years    |   14930126 (27.3%)   |  14930126 (27.3%)  |
#   |-  80+ years      |   11105538 (20.3%)   |  11105538 (20.3%)  |


##### 4.1.3 种族-race #####
table(paper.data$RIDRETH3)
#recode_factor 类似编码which 因子转换
race <- recode_factor(paper.data$RIDRETH3, 
                      `1` = 'Mexican American',
                      `2` = 'Other Hispanic',
                      `3` = 'Non-Hispanic White',
                      `4` = 'Non-Hispanic Black',
                      `6` = 'Other/multiracial',
                      `7` = 'Other/multiracial'
)
paper.data$race <- race
table(paper.data$race)

tab1 <- tableby(Type ~ factor(race),
                data = paper.data, weights=WTDRD1,
                subset = WTDRD1>0) 

summary(tab1, text=TRUE)

# |                      | include (N=54660094) | Total (N=54660094) |
#   |:---------------------|:--------------------:|:------------------:|
#   |factor(race)          |                      |                    |
#   |-  Mexican American   |    1921207 (3.5%)    |   1921207 (3.5%)   |
#   |-  Other Hispanic     |    2088012 (3.8%)    |   2088012 (3.8%)   |
#   |-  Non-Hispanic White |   43095361 (78.8%)   |  43095361 (78.8%)  |
#   |-  Non-Hispanic Black |    4665975 (8.5%)    |   4665975 (8.5%)   |
#   |-  Other/multiracial  |    2889538 (5.3%)    |   2889538 (5.3%)   |
##### 4.1.4 复杂转换-饮酒 alq.group #####
# 仅衍生有无饮酒
# 区分男性和女性不同的饮酒标准，drinks per day (dkspd) ：
# 衍生方式：Paper01-美国成年人膳食中类胡萝素与认知功能的关系, NHANES 2011-2014
# 衍生结果为：Non-drinker, 1-5 drinks/month, 5-10 drinks/month, 10+ drinks/month
# ALQ120Q, 具体的数值;  777:refused, 999:don't know
# ALQ120U, 1:week, 2:month, 3:year, 7:refused, 9:don't know

# aql110 一生当中都至少喝过十二次酒
# alq101  过去十二个月 至少十二次酒

# 使用 plyr 包的 ddply 快速看一下数据 可以快速查看表格数据进行统计出来
# 根据数据的alq101 和alq110 进行分别对应去取数量 根据seqn 列出不同条件下统计的的数量
ddply(paper.data, .(ALQ101, ALQ110), summarize, n = length(SEQN))
ddply(alq_g, .(ALQ101, ALQ110), summarize, n = length(SEQN))
####
# 统一单位为月
# week-month: *4; year-month:/12
ori.alq.unit <- paper.data$ALQ120U
trans.unit.month <- ifelse(ori.alq.unit==1,4,
                           ifelse(ori.alq.unit==3,1/12,
                                  ifelse(ori.alq.unit==7| ori.alq.unit==9,NA,1)))

paper.data$trans.unit.month <- trans.unit.month


# 统一数值为月饮酒量-quantity
ori.alq.quantity <- paper.data$ALQ120Q
trans.quantity.month <- ifelse(ori.alq.quantity >= 0,  
                               ori.alq.quantity * trans.unit.month, NA)
paper.data$trans.quantity.month <- trans.quantity.month
#View(paper.data[,c('SEQN','ALQ101','ALQ110', 'ALQ120Q', 'ALQ120U', 'trans.quantity.month', 'trans.unit.month')])

# 将饮酒量的结果划分为4档：Non-drinker, 1-5 drinks/month, 5-10 drinks/month, 10+ drinks/month
# 首先根据转化得到的trans.quantity.month，将饮酒量划分3档
# 利用 paper.data$ALQ101
alq101 <- paper.data$ALQ101
trans.quantity.month.factor <- ifelse(trans.quantity.month>=1 & trans.quantity.month<5,'1-5 drinks/month',
                                      ifelse(trans.quantity.month >=5 & trans.quantity.month<10,'5-10 drinks/month',
                                      ifelse(trans.quantity.month>=10,'10+ drinks/month','wait')))

paper.data$trans.quantity.month.factor <- trans.quantity.month.factor
ddply(paper.data, .(ALQ101, trans.quantity.month.factor), summarise, n = length(SEQN))
# 单位转换需要时刻注意 需要多使用view进行查看数据 与官网进行指控 看了视频后这段与这篇论文生成的数据不一致 应该是作者计算没有进行单位转换,

# alq101 是1 但是trans.quantity.month.factor 没有值 NA 和wait 直接改掉
alq101Nawait <- which((trans.quantity.month.factor=="wait"| is.na(trans.quantity.month.factor))& alq101==1 )
table(alq101Nawait)
trans.quantity.month.factor[alq101Nawait] <-'1-5 drinks/month'
table(trans.quantity.month.factor)

paper.data$trans.quantity.month.factor <- trans.quantity.month.factor
ddply(paper.data, .(ALQ101, trans.quantity.month.factor), summarise, n = length(SEQN))

index.less.1 <- which(alq101 == 2)
trans.quantity.month.factor[index.less.1] <- 'Non-drinker'
table(trans.quantity.month.factor)
paper.data$alq.group <- trans.quantity.month.factor
table(paper.data$alq.group)

tab1 <- tableby(Type ~ factor(alq.group) ,
                data = paper.data, weights=WTDRD1,
                subset = WTDRD1>0) 

summary(tab1, text=TRUE)

# |                     | include (N=54660094) | Total (N=54660094) |
#   |:--------------------|:--------------------:|:------------------:|
#   |factor(alq.group)    |                      |                    |
#   |-  N-Miss            |        286384        |       286384       |
#   |-  1-5 drinks/month  |   26113325 (48.0%)   |  26113325 (48.0%)  |
#   |-  10+ drinks/month  |   10742749 (19.8%)   |  10742749 (19.8%)  |
#   |-  5-10 drinks/month |    2778425 (5.1%)    |   2778425 (5.1%)   |
#   |-  Non-drinker       |   14739211 (27.1%)   |  14739211 (27.1%)  |
#   |factor(RIAGENDR)     |                      |                    |
#   |-  1                 |   25204539 (46.1%)   |  25204539 (46.1%)  |
#   |-  2                 |   29455555 (53.9%)   |  29455555 (53.9%)  |

tab1 <- tableby(Type ~ factor(Sex) + factor(age.group)+ factor(race)+ factor(alq.group),
                data = paper.data, weights=WTDRD1,
                subset = WTDRD1>0) 

summary(tab1, text=TRUE)
# |                      | include (N=54660094) | Total (N=54660094) |
#   |:---------------------|:--------------------:|:------------------:|
#   |factor(Sex)           |                      |                    |
#   |-  female             |   29455555 (53.9%)   |  29455555 (53.9%)  |
#   |-  male               |   25204539 (46.1%)   |  25204539 (46.1%)  |
#   |factor(age.group)     |                      |                    |
#   |-  60-69 years        |   28624430 (52.4%)   |  28624430 (52.4%)  |
#   |-  70-79 years        |   14930126 (27.3%)   |  14930126 (27.3%)  |
#   |-  80+ years          |   11105538 (20.3%)   |  11105538 (20.3%)  |
#   |factor(race)          |                      |                    |
#   |-  Mexican American   |    1921207 (3.5%)    |   1921207 (3.5%)   |
#   |-  Other Hispanic     |    2088012 (3.8%)    |   2088012 (3.8%)   |
#   |-  Non-Hispanic White |   43095361 (78.8%)   |  43095361 (78.8%)  |
#   |-  Non-Hispanic Black |    4665975 (8.5%)    |   4665975 (8.5%)   |
#   |-  Other/multiracial  |    2889538 (5.3%)    |   2889538 (5.3%)   |
#   |factor(alq.group)     |                      |                    |
#   |-  N-Miss             |        286384        |       286384       |
#   |-  1-5 drinks/month   |   26113325 (48.0%)   |  26113325 (48.0%)  |
#   |-  10+ drinks/month   |   10742749 (19.8%)   |  10742749 (19.8%)  |
#   |-  5-10 drinks/month  |    2778425 (5.1%)    |   2778425 (5.1%)   |
#   |-  Non-drinker        |   14739211 (27.1%)   |  14739211 (27.1%)  |

##### 4.1.5 Total 卡路里-total.calories #####
# day1、day2 都有，则取2天平均，否则，取第1天的值
#apply 函数可以对行和列进行批量操作 前面给出两列 如果2列有值则取平均 否则就只取第一列 1是对每一行去取平均值
#View(paper.data[,c('DR1TKCAL', 'DR2TKCAL')])
# 先计算2天的平均值，其中第二天为 NA 的，平均值也是 NA
total.calories <- apply(paper.data[,c('DR1TKCAL', 'DR2TKCAL')], 1, mean)
paper.data$total.calories <- total.calories
#View(paper.data[,c('DR1TKCAL', 'DR2TKCAL', 'total.calories')])

#处理掉na 这是拿到表的行数 直接赋值
day.2.na.index<-which(is.na(paper.data$DR2TKCAL))
# 这个是取的data.frame的 列的行值集合
total.calories[day.2.na.index]<-paper.data$DR1TKCAL[day.2.na.index]

paper.data$total.calories <- total.calories

##### 4.1.6 Total 饮食中 L&Z-total.dr.LZ#####
# https://wwwn.cdc.gov/Nchs/Nhanes/2011-2012/DR2TOT_G.htm#DR2TLZ
# day1、day2 都有，则取2天平均，否则，取第一天的值
# 先计算2天的平均值，其中第二天为 NA 的，平均值也是 NA
# apply apply(X, MARGIN, FUN, ...)
# X: 数组、矩阵、数据框，数据至少是二维的
# MARGIN: 按行计算或按列计算，1 表示按行，2 表示按列
# FUN: 自定义的调用函数

# 下面是就是两列n 行 1代表行 就是取每行平均值返回 

total.dr.LZ <- apply(paper.data[,c('DR1TLZ', 'DR2TLZ')], 1, mean)
paper.data$total.dr.LZ <- total.dr.LZ
#View(paper.data[,c('DR1TLZ', 'DR2TLZ', 'total.dr.LZ')])

# 把第二天为 NA 的值计算出来，用第一天的值作为平均值
day.2.na.index <- which(is.na(paper.data$total.dr.LZ))
total.dr.LZ[day.2.na.index] <- paper.data$DR1TLZ[day.2.na.index]
total.dr.LZ[day.2.na.index]

# 添加到paper.data中，再次确认
paper.data$total.dr.LZ <- total.dr.LZ
# View(paper.data[,c('DR1TLZ', 'DR2TLZ', 'total.dr.LZ')])


##### 4.1.7 认知评分 #####
#CERAD 词汇第1次即刻记忆测验	CFDCST1
#CERAD 词汇第2次即刻记忆测验	CFDCST2
#CERAD 词汇第3次即刻记忆测验	CFDCST3
#CERAD 3次词汇即刻记忆测验总分	CERAD_1+CERAD_2+CERAD_3

# colnames(paper.data)
CERAD.total <- apply(paper.data[,c('CFDCST1', 'CFDCST2', 'CFDCST3')], 1, sum)
paper.data$CERAD.total <- CERAD.total
#View(paper.data[,c('CFDCST1', 'CFDCST2', 'CFDCST3','CERAD.total')])











